{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Stochastic examples\n",
    "\n",
    "\n",
    "This example is designed to show how to use the stochatic optimization\n",
    "algorithms for descrete and semicontinous measures from the POT library.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Kilian Fatras <kilian.fatras@gmail.com>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import matplotlib.pylab as pl\n",
    "import numpy as np\n",
    "import ot\n",
    "import ot.plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "COMPUTE TRANSPORTATION MATRIX FOR SEMI-DUAL PROBLEM\n",
    "############################################################################\n",
    "############################################################################\n",
    " DISCRETE CASE:\n",
    "\n",
    " Sample two discrete measures for the discrete case\n",
    " ---------------------------------------------\n",
    "\n",
    " Define 2 discrete measures a and b, the points where are defined the source\n",
    " and the target measures and finally the cost matrix c.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_source = 7\n",
    "n_target = 4\n",
    "reg = 1\n",
    "numItermax = 1000\n",
    "\n",
    "a = ot.utils.unif(n_source)\n",
    "b = ot.utils.unif(n_target)\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "X_source = rng.randn(n_source, 2)\n",
    "Y_target = rng.randn(n_target, 2)\n",
    "M = ot.dist(X_source, Y_target)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call the \"SAG\" method to find the transportation matrix in the discrete case\n",
    "---------------------------------------------\n",
    "\n",
    "Define the method \"SAG\", call ot.solve_semi_dual_entropic and plot the\n",
    "results.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.55553509e-02 9.96395660e-02 1.76579142e-02 4.31178196e-06]\n",
      " [1.21640234e-01 1.25357448e-02 1.30225078e-03 7.37891338e-03]\n",
      " [3.56123975e-03 7.61451746e-02 6.31505947e-02 1.33831456e-07]\n",
      " [2.61515202e-02 3.34246014e-02 8.28734709e-02 4.07550428e-04]\n",
      " [9.85500870e-03 7.52288517e-04 1.08262628e-02 1.21423583e-01]\n",
      " [2.16904253e-02 9.03825797e-04 1.87178503e-03 1.18391107e-01]\n",
      " [4.15462212e-02 2.65987989e-02 7.23177216e-02 2.39440107e-03]]\n"
     ]
    }
   ],
   "source": [
    "method = \"SAG\"\n",
    "sag_pi = ot.stochastic.solve_semi_dual_entropic(a, b, M, reg, method,\n",
    "                                                numItermax)\n",
    "print(sag_pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SEMICONTINOUS CASE:\n",
    "\n",
    "Sample one general measure a, one discrete measures b for the semicontinous\n",
    "case\n",
    "---------------------------------------------\n",
    "\n",
    "Define one general measure a, one discrete measures b, the points where\n",
    "are defined the source and the target measures and finally the cost matrix c.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_source = 7\n",
    "n_target = 4\n",
    "reg = 1\n",
    "numItermax = 1000\n",
    "log = True\n",
    "\n",
    "a = ot.utils.unif(n_source)\n",
    "b = ot.utils.unif(n_target)\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "X_source = rng.randn(n_source, 2)\n",
    "Y_target = rng.randn(n_target, 2)\n",
    "M = ot.dist(X_source, Y_target)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call the \"ASGD\" method to find the transportation matrix in the semicontinous\n",
    "case\n",
    "---------------------------------------------\n",
    "\n",
    "Define the method \"ASGD\", call ot.solve_semi_dual_entropic and plot the\n",
    "results.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.88833283 7.64041833 3.93000933 2.68489048 1.42837354 3.25840738\n",
      " 2.80033951] [-2.50038759 -2.4083026  -0.96389053  5.87258072]\n",
      "[[2.49326139e-02 1.01118047e-01 1.68018025e-02 4.67918477e-06]\n",
      " [1.20543018e-01 1.29218840e-02 1.25860644e-03 8.13363473e-03]\n",
      " [3.52425849e-03 7.83826265e-02 6.09501106e-02 1.47316769e-07]\n",
      " [2.62727985e-02 3.49290291e-02 8.11998888e-02 4.55426386e-04]\n",
      " [9.00986942e-03 7.15412954e-04 9.65318348e-03 1.23478677e-01]\n",
      " [1.98446848e-02 8.60145164e-04 1.67017745e-03 1.20482135e-01]\n",
      " [4.16774129e-02 2.77550575e-02 7.07529364e-02 2.67173611e-03]]\n"
     ]
    }
   ],
   "source": [
    "method = \"ASGD\"\n",
    "asgd_pi, log_asgd = ot.stochastic.solve_semi_dual_entropic(a, b, M, reg, method,\n",
    "                                                           numItermax, log=log)\n",
    "print(log_asgd['alpha'], log_asgd['beta'])\n",
    "print(asgd_pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the results with the Sinkhorn algorithm\n",
    "---------------------------------------------\n",
    "\n",
    "Call the Sinkhorn algorithm from POT\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.55535622e-02 9.96413843e-02 1.76578860e-02 4.31043335e-06]\n",
      " [1.21640742e-01 1.25369034e-02 1.30234529e-03 7.37715259e-03]\n",
      " [3.56096458e-03 7.61460101e-02 6.31500344e-02 1.33788624e-07]\n",
      " [2.61499607e-02 3.34255577e-02 8.28741973e-02 4.07427179e-04]\n",
      " [9.85698720e-03 7.52505948e-04 1.08291770e-02 1.21418473e-01]\n",
      " [2.16947591e-02 9.04086158e-04 1.87228707e-03 1.18386011e-01]\n",
      " [4.15442692e-02 2.65998963e-02 7.23192701e-02 2.39370724e-03]]\n"
     ]
    }
   ],
   "source": [
    "sinkhorn_pi = ot.sinkhorn(a, b, M, reg)\n",
    "print(sinkhorn_pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PLOT TRANSPORTATION MATRIX\n",
    "#############################################################################\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot SAG results\n",
    "----------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAExZJREFUeJzt3X+wpQV93/H3hwUVAhHj3ihhxbVqd4pawdxgGqwa/IXE/JqYalJ/RdutrVhpTa0mmY7WaZImUyWd2qRbNcaIEo06k+aHhQYYytQfvasbhh8ygw66rCAXCAoUsCzf/vGcbe/c7u49u3vO+e6e837NnNl773nOeb7nwH3f5zznOeekqpAkzd5x3QNI0qIywJLUxABLUhMDLElNDLAkNTHAktTEAGvuJPn7SS47hOXfkOSaCa37liQvnsR1HYuSnJHkviSbumc5FhhgzZ2quqSqXto9x+FI8ookX0pyf5K7klySZMvovF8Zxe2+JA8m2bvm++tnMNuGf1yq6ptVdXJV7T2M639GksuS3J3kniQ7k1ywbpmnJHkkye/u5/JJcmGSa5P8ryS3J7kqyasPdZZZMcDSUSLJK4GPAxcDm4FnAA8B1yR5XFX9+ihuJwNvBj6/7/uqekbf5IMkxx/hVfwX4HLgicAPAv8U+O66ZV4H/DXwqiSPXnfevwcuAt4OPB44Hfg14PwjnGt6qsqTp5mdgH8J7AHuBW4CXjT6+XHAO4GvAXcBnwR+YHTeVqCAXwJ2M/wCvhn4EeBa4B7gP6xZxxuAaw4yw+OBP2H45f4S8N59y69Z1/Frlr8K+Aejr58KXDGa8U7gEuDUNcveArz4MO6XAN8A3rHu58cB1wH/et3PD3obD/N+O+BtA/4QeAR4ALgPeMea638T8E3g6rX3H/ADwK3AT46u42TgZuB1+5l18+hyp25wH30N+MfAt4FXrjnvbwJ7geXu/8cP5eQWsGYmyTbgQuBHquoU4GUMwQJ4K/AzwAuAH2KIxQfWXcVzgacDr2LYSvxV4MUMW4p/L8kLxhzlA8CDwGnAG0ensW8G8BujGf8W8CTg3WNdMPnFJNce4OxtwBnAp9b+sKoeAT4NvOQQZlxv3PvtgLetql7LENmfrGGL+7fWXP8LRsu/bN3sdzPct/85yQ8C7wd2VdVH9zPjXQxx/liSn0nyhP0s8zxgC3Apwx/o16857zxgd1WtbHhvHEUMsGZpL/Bo4MwkJ1TVLVX1tdF5bwZ+tapuraqHGH7xX7nuYe17q+rBqroMuB/4RFXdUVV7gP8OnL3RAKMnh34O+FdVdX9VXQf8wbg3oKpurqrLq+qhqloF3scQoHEu+/Gq+tsHOHvz6N/b9nPebWvOPxxj3W9HcNvePbovH1h/xmidnwL+ErgA+Ef7u4IaNmN/nOEP8r8DbktydZKnr1ns9cBfVNVfM+yqOX8Udhjun9vXXmeSW0f7kh9M8uQxbsfMGWDNTFXdzLCP7t3AHUkuTfJDo7OfDHx29AtzD3AjQ7DXbgl9e83XD+zn+5PXr3PdE1e/BywxPDzevWaxb4x7G5I8YTT3niTfBT7GkcVxnztH/562n/NOW3P+4RjrfjuC27Z7g/N3AM8EPlJVdx1oodEf3wur6qkM/z/cD3x0NNuJwM8z7Bahqj7PsEX+i6OL38W6+66qtozmfzTD1v1RxwBrpkZbgc9j+AUr4N+OztoNvLyqTl1zesxoK+1I1vd/n7iqqjcDq8DDDA+v9zljzdf3j/49ac3Pnrjm618fzf2sqvp+4DVM5pf7Job9pT+/9odJjmPYYv/LCaxjIxvdtgO9deIB31Jx9IhjB0NI/0mSp40zSFXtZthV9MzRj34W+H7gP46Obrid4Um2fbshrgC2JFke5/qPFgZYM5NkW5LzRs9eP8iw9fXI6OzfA/7NvoeKSZaS/PSkZ6jh8KjPAO9OclKSM1mzL3H00HsP8Jokm5K8keHJqX1OYXgS6jtJTgf+xYTmKuCXgV8b7St+TJInAh9kCM/7J7GeDWx0274N/I1DvM5fYQj0G4HfBj66v2OEkzwuyXuSPC3JcUk2jy7zhdEirwc+DDwLOGt0Ohd4dpJnVdVNwH8CLk3ykiQnjtbzY4c470wZYM3So4HfZHg4fTvDoUbvGp33OwxHJlyW5F6GX7znTmmOCxkedt8OfAT4/XXn/0OG+NzF8ETV/1hz3nuA5wDfAf6MIeZjGb1A5IDH61bVHwGvBf7ZaN03ACcC5x7sofsEbXTbfoPhD8Q9SX55oytL8sPAP2c46mEvw6OdYjjaZb3vMRxB8d8Yjk65juEQvDeM/hi8CLi4qm5fc9oJfI7/9wf0LQyHor0PuJvhEcV7GZ58/OZY98CMZXQIhyRpxtwClqQmBliSmhhgSWpigCWpyZG+eYaOcZs3b66tW7d2jyHNlZ07d95ZVUsbLWeAF9zWrVtZWTmmXj4vHfWSjPXqSndBSFITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNju8eQM1uugle+MLuKbTIzjoLLr64e4oWbgFLUhO3gBfdtm1w1VXdU0gLyS1gSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJqmq7hnUKMm9wE3dc0zYZuDO7iGmwNt17NhWVadstNDxs5hER7Wbqmq5e4hJSrIyb7cJvF3HkiQr4yznLghJamKAJamJAdaO7gGmYB5vE3i7jiVj3SafhJOkJm4BS1ITAyxJTQzwgkpyfpKbktyc5J3d80xCkg8nuSPJdd2zTEqSJyW5MskNSa5P8rbumSYhyWOSfCnJX41u13u6Z5qUJJuSfCXJn260rAFeQEk2AR8AXg6cCfxCkjN7p5qIjwDndw8xYQ8Db6+qM4EfBd4yJ/+tHgLOq6pnA2cB5yf50eaZJuVtwI3jLGiAF9M5wM1V9fWq+h5wKfDTzTMdsaq6Gri7e45JqqrbqurLo6/vZfjFPr13qiNXg/tG354wOh3zRwQk2QL8BPDBcZY3wIvpdGD3mu9vZQ5+qeddkq3A2cAXeyeZjNFD9V3AHcDlVTUPt+ti4B3AI+MsbIClY0CSk4FPAxdV1Xe755mEqtpbVWcBW4Bzkjyze6YjkeQVwB1VtXPcyxjgxbQHeNKa77eMfqajUJITGOJ7SVV9pnueSauqe4ArOfb3358L/FSSWxh2652X5GMHu4ABXkz/E3h6kqckeRTwauBPmmfSfiQJ8CHgxqp6X/c8k5JkKcmpo69PBF4CfLV3qiNTVe+qqi1VtZXhd+qKqnrNwS5jgBdQVT0MXAj8V4YndT5ZVdf3TnXkknwC+DywLcmtSd7UPdMEnAu8lmFratfodEH3UBNwGnBlkmsZNggur6oND9uaN74UWZKauAUsSU2m8obsmzdvrq1bt07jqjVhO3fuvLOqlrrnOFIvfOlvHtZDuZe9/+pJj3JQV77unJmuD6C+Mtu9S5c/8qnMdIXHsKkEeOvWraysjPWG8GqW5BvdM0iLyl0QktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUZKwAz+MHOEpStw0DPMcf4ChJrcbZAp7LD3A8FBddNJwkaZLGeTOe/X2A43PXL5RkO7Ad4IwzzpjIcEeLXbu6J5A0jyb2JFxV7aiq5apaXlo65t/dUJKmbpwA+wGOkjQF4wTYD3CUpCnYcB9wVT2cZN8HOG4CPjwPH+AoSd3G+kSMqvpz4M+nPIskLRRfCSdJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSk7FeiCEd7a746IcO63IXPP9nJzzJwdXXvzrT9QFs8s2xjlpuAUtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNNgxwkg8nuSPJdbMYSJIWxThbwB8Bzp/yHJK0cDYMcFVdDdw9g1kkaaG4D1iSmkwswEm2J1lJsrK6ujqpq5WkuTWxAFfVjqparqrlJd9/VJI25C4ISWoyzmFonwA+D2xLcmuSN01/LEmafxt+JFFV/cIsBpGkReMuCElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJarLhK+GkY8HLn/Zjh3W5b/7hSROe5OAe+NbyTNcH8PS3fnHm69R43AKWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmozzqchPSnJlkhuSXJ/kbbMYTJLm3TjvBfEw8Paq+nKSU4CdSS6vqhumPJskzbUNt4Cr6raq+vLo63uBG4HTpz2YJM27Q9oHnGQrcDbw/729UpLtSVaSrKyurk5mOkmaY2MHOMnJwKeBi6rqu+vPr6odVbVcVctLS0uTnFGS5tJYAU5yAkN8L6mqz0x3JElaDOMcBRHgQ8CNVfW+6Y8kSYthnC3gc4HXAucl2TU6XTDluSRp7m14GFpVXQNkBrNI0kLxlXCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNRnn/YClo96Df/fMw7rcYz8521+Bx7/x2zNdn45ubgFLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTcb5VOTHJPlSkr9Kcn2S98xiMEmad+O8EP4h4Lyqui/JCcA1Sf6iqr4w5dkkaa6N86nIBdw3+vaE0ammOZQkLYKx9gEn2ZRkF3AHcHlVfXE/y2xPspJkZXV1ddJzStLcGSvAVbW3qs4CtgDnJHnmfpbZUVXLVbW8tLQ06Tklae4c0lEQVXUPcCVw/nTGkaTFMc5REEtJTh19fSLwEuCr0x5MkubdOEdBnAb8QZJNDMH+ZFX96XTHkqT5N85RENcCZ89gFklaKL4STpKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQm47wSTjrqnXT9bYd1uUft+daEJzm447+wZabrA/izb+2a+To1HreAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCZjBzjJpiRfSeIHckrSBBzKFvDbgBunNYgkLZqxApxkC/ATwAenO44kLY5xt4AvBt4BPHKgBZJsT7KSZGV1dXUiw0nSPNswwEleAdxRVTsPtlxV7aiq5apaXlpamtiAkjSvxtkCPhf4qSS3AJcC5yX52FSnkqQFsGGAq+pdVbWlqrYCrwauqKrXTH0ySZpzHgcsSU0O6SOJquoq4KqpTCJJC8YtYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaHNILMaSj1f9+8uG9AVT2fGvCkxzc3j23zXR9AN955IGZru9xM13bsc0tYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJanJWC9FHn0k/b3AXuDhqlqe5lCStAgO5b0gfryq7pzaJJK0YNwFIUlNxg1wAZcl2Zlk+/4WSLI9yUqSldXV1clNKElzatwAP6+qngO8HHhLkuevX6CqdlTVclUtLy0d3lsDStIiGSvAVbVn9O8dwGeBc6Y5lCQtgg0DnOT7kpyy72vgpcB10x5MkubdOEdBPAH4bJJ9y3+8qj431akkaQFsGOCq+jrw7BnMIkkLxcPQJKmJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpyaG8Ibt01LrzWSce1uUee/IPT3iSg9v9+odnuj6AVz1100zXd9kDM13dMc0tYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJajJWgJOcmuSPk3w1yY1J/s60B5OkeTfuS5F/B/hcVb0yyaOAk6Y4kyQthA0DnOSxwPOBNwBU1feA7013LEmaf+PsgngKsAr8fpKvJPlgku+b8lySNPfGCfDxwHOA362qs4H7gXeuXyjJ9iQrSVZWV1cnPGavs84aTpI0SePsA74VuLWqvjj6/o/ZT4CragewA2B5ebkmNuFR4OKLuyeQNI823AKuqtuB3Um2jX70IuCGqU4lSQtg3KMg3gpcMjoC4uvAL01vJElaDGMFuKp2ActTnkWSFoqvhJOkJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCapmvz75iRZBb4x8SvWNDy5qpa6h5AW0VQCLEnamLsgJKmJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCb/B6HXs8MRx/3SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, sag_pi, 'semi-dual : OT matrix SAG')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot ASGD results\n",
    "-----------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE3tJREFUeJzt3X+wpQV93/H3xwUFXRKiezXIgmus2YbYBPSKSbFRUQwqMZmJrfiL+KPdOhUHWlOLSdrRZvKj7YzFjE7SrTHEilKMOskkJoGJMIapP3JXN4Qf0hIGZYnARYr8COKw++0fz9nMze3u3rO759zv7jnv18yZvfc+zznP95xl3zz3Oc85J1WFJGn9Pa57AEmaVwZYkpoYYElqYoAlqYkBlqQmBliSmhhgHdWSvCHJVQex/puTXDehbd+e5GWTuK2jUZJTkzyUZEP3LEcrA6yjWlVdXlUv757jUCQ5L8mXkzyc5FtJLk+yebTsF0ZxeyjJd5LsXvH9jesw25r/c6mqb1TVxqrafRjbuSzJY0lOWvXzE5N8JMldSR5M8r+TXLJieZJcmOT6JH87Wu/aJOevWOfa0WP3YJIHkuxIckmSJxzqvJNmgKUGSV4DfBy4FNgE/DDwKHBdku+rql8dxW0j8HbgC3u/r6of7pt8kOSYCdzGk4CfBb4NvHHV4v8KbAR+CPhe4NXArSuW/wZwMfAu4CnAycAvAeeuup0Lq+oE4KTRuucDn02Sw51/IqrKi5epXIB/B9wJPAjcArx09PPHAZcAfw18C7gSePJo2RaggLcAdwD/lyFAzweuB+4HPrhiG28GrjvADE8B/gB4APgy8Mt711+xrWNWrH8t8M9HXz8L+NxoxnuBy4ETV6x7O/CyQ3hcAnwdePeqnz8OuAH4j6t+fsD7eIiP237vG/A/gD3AI8BDwLtX3P7bgG8An1/5+AFPBnYBPzW6jY0MwbzgADNfMJr1IuCGVctuAH5mP9f7QWA3sLjGY/J3f5crfnYq8LfAed3/PqrKPWBNR5KtwIXA82vYA/lJhmABvBP4GeBFwNMZYvGhVTfxAuDZwGsZ9hJ/EXgZw57iP0vyojFH+RDwHYY9oLeOLmPfDeDXRjP+EHAK8N6xrpi8Psn1+1m8lSEEn1z5w6raA3wKOOcgZlxt3Mdtv/etqt7EENmfqmGP+z+vuP0Xjdb/yVWz38fw2P73JE9l2IPdWVUfPcCsPwd8ArgC+IdJnrdi2ReBX0nyliTPXnW9s4E7qmppjcfi/1NV3wCWgH9ysNedBgOsadkNPAE4LcmxVXV7Vf31aNnbgV+sql1V9SjDP/zXrPq19per6jtVdRXwMPCJqrqnqu4E/hw4Y60BRk8O/SzwH6rq4aq6Afjdce9AVd1aVVdX1aNVtQy8nyFA41z341X1I/tZvGn05zf3seybK5YfirEet8O4b+8dPZaPrF4w2uYngT8DXgn8y/3dSJJTgZcAH6+qu0fXuWDFKu9k2Cu/ELgpya1JXjFatgm4a9Xt7Upy/+iY7zPWuA9/w7DH3s4Aayqq6laGY3TvBe5JckWSp48WPwP4zOgfzP3AzQzBftqKm7h7xdeP7OP7jau3ueqJq98CFhh+Pb5jxWpfH/c+JHnaaO47kzwAfIzDi+Ne947+PGkfy05asfxQjPW4HcZ9u2ON5duB5wCXVdW3DrDem4Cbq2rn6PvLgdcnORagqh6p4Tj48xgOI10JfDLJkxkOm/y9x66qNo/mfwLD3v2BnAzct8Y668IAa2pGe4EvZAhuAf9ptOgO4BVVdeKKy3GjvbTD2d7fPXFVVW8HloHHGH693uvUFV8/PPrziSt+9v0rvv7V0dz/qKq+h+GJokk8eXMLw/HSf7ryh0kex7DH/mcT2MZa1rpv+3ubxP2+feLoN47twEeBf5XkHxxg+xcAPzA6e+Euhj3wTQx7zn9/g1UPjOZ9EvBMhmPXm5MsHuD29zfjKcDzGH4baGeANRVJtiY5e3TKz3cY9r72jBb/FsPxvWeM1l1I8tOTnqGG06M+Dbw3yROTnMZw3HHv8mWGJwnfmGRDkrcyPDm11wkMT0J9O8nJwL+d0FwF/DzwS6Njxccl+X7gw8D3MBw/nba17tvdwA8c5G3+AkOg3wr8F+Cj+zpHOMmPMzzOZwKnjy7PYTgr5ILROv8+yfOTPD7JcQxP1N0P3FJVtwD/DbgiyTlJjh9t5x/vb7DR3/+LgN9neDL2swd536bCAGtangD8OsOv03cBTwXeM1r2AYYzE65K8iDDEy4vmNIcFzL82n0XcBnwO6uW/wuG+HyL4Ymq/7Vi2fuA5zKcJvVHDDEfS4YXiOz3fN2q+p8Mv4b/69G2bwKOB85a41f3SVnrvv0aw/8g7k/y82vd2OgJtH/DcNbDbobfdorhbJfVfg74/ar6q6q6a++F4b+L80aHGYrh7+pehmO25wCvqqqHRrfxDoZT0d7PcDhhF8MZLq9leAJxrw+O/hu7m+FJyU8B546e8GyX0akZkqR15h6wJDUxwJLUxABLUhMDLElNDvsNNXR027RpU23ZsqV7DGmm7Nix496qWlhrPQM857Zs2cLS0kG/pF7SASQZ6xWXHoKQpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoc0z2Amt1yC7z4xd1TaJ6dfjpcemn3FC3cA5akJu4Bz7utW+Haa7unkOaSe8CS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNUlXdM6hRkgeBW7rnmLBNwL3dQ0yB9+vosbWqTlhrpWPWYxId0W6pqsXuISYpydKs3Sfwfh1NkiyNs56HICSpiQGWpCYGWNu7B5iCWbxP4P06mox1n3wSTpKauAcsSU0MsCQ1McBzKsm5SW5JcmuSS7rnmYQkH0lyT5IbumeZlCSnJLkmyU1JbkxyUfdMk5DkuCRfTvKXo/v1vu6ZJiXJhiRfTfKHa61rgOdQkg3Ah4BXAKcBr0tyWu9UE3EZcG73EBP2GPCuqjoN+DHgHTPyd/UocHZV/ShwOnBukh9rnmlSLgJuHmdFAzyfzgRurarbquq7wBXATzfPdNiq6vPAfd1zTFJVfbOqvjL6+kGGf9gn9051+Grw0OjbY0eXo/6MgCSbgVcBHx5nfQM8n04G7ljx/S5m4B/1rEuyBTgD+FLvJJMx+lV9J3APcHVVzcL9uhR4N7BnnJUNsHQUSLIR+BRwcVU90D3PJFTV7qo6HdgMnJnkOd0zHY4k5wH3VNWOca9jgOfTncApK77fPPqZjkBJjmWI7+VV9enueSatqu4HruHoP35/FvDqJLczHNY7O8nHDnQFAzyf/gJ4dpJnJnk8cD7wB80zaR+SBPht4Oaqen/3PJOSZCHJiaOvjwfOAb7WO9Xhqar3VNXmqtrC8G/qc1X1xgNdxwDPoap6DLgQ+FOGJ3WurKobe6c6fEk+AXwB2JpkV5K3dc80AWcBb2LYm9o5uryye6gJOAm4Jsn1DDsEV1fVmqdtzRpfiixJTdwDlqQmU3lD9k2bNtWWLVumcdOasB07dtxbVQvdcxyul7zs1w/pV7mXf+Dzkx7lgK59w/PWdXsAe65f30OrV+/5ZNZ1g0exqQR4y5YtLC2N9Ybwapbk690zSPPKQxCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktRkrADP4gc4SlK3NQM8wx/gKEmtxtkDnskPcDwYF188XCRpksZ5M559fYDjC1avlGQbsA3g1FNPnchwR4qdO7snkDSLJvYkXFVtr6rFqlpcWDjq391QkqZunAD7AY6SNAXjBNgPcJSkKVjzGHBVPZZk7wc4bgA+Mgsf4ChJ3cb6RIyq+izw2SnPIklzxVfCSVITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktRkrBdiSEe6T3/0g4d0vde9+PUTnuTA9tz2f9Z1ewAbnvbUdd+mxuMesCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktRkzQAn+UiSe5LcsB4DSdK8GGcP+DLg3CnPIUlzZ80AV9XngfvWYRZJmiseA5akJhMLcJJtSZaSLC0vL0/qZiVpZk0swFW1vaoWq2pxYWFhUjcrSTPLQxCS1GSc09A+AXwB2JpkV5K3TX8sSZp9a34kUVW9bj0GkaR54yEISWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqsuYr4aSjwfk/+NJDut6uy4+f8CQH9shtz1/X7QE8611fXPdtajzuAUtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNxvlU5FOSXJPkpiQ3JrloPQaTpFk3zntBPAa8q6q+kuQEYEeSq6vqpinPJkkzbc094Kr6ZlV9ZfT1g8DNwMnTHkySZt1BHQNOsgU4A/jSPpZtS7KUZGl5eXky00nSDBs7wEk2Ap8CLq6qB1Yvr6rtVbVYVYsLCwuTnFGSZtJYAU5yLEN8L6+qT093JEmaD+OcBRHgt4Gbq+r90x9JkubDOHvAZwFvAs5OsnN0eeWU55KkmbfmaWhVdR2QdZhFkuaKr4STpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqck47wcsHfEefeFph3S9jVeu7z+B495w37puT0c294AlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJuN8KvJxSb6c5C+T3JjkfesxmCTNunFeCP8ocHZVPZTkWOC6JH9cVV+c8mySNNPG+VTkAh4afXvs6FLTHEqS5sFYx4CTbEiyE7gHuLqqvrSPdbYlWUqytLy8POk5JWnmjBXgqtpdVacDm4EzkzxnH+tsr6rFqlpcWFiY9JySNHMO6iyIqrofuAY4dzrjSNL8GOcsiIUkJ46+Ph44B/jatAeTpFk3zlkQJwG/m2QDQ7CvrKo/nO5YkjT7xjkL4nrgjHWYRZLmiq+Ek6QmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJanJOK+Ek454x//VrkO63rF33T3hSQ7smD9/+rpuD+CP/mbnum9T43EPWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWoydoCTbEjy1SR+IKckTcDB7AFfBNw8rUEkad6MFeAkm4FXAR+e7jiSND/G3QO+FHg3sGd/KyTZlmQpydLy8vJEhpOkWbZmgJOcB9xTVTsOtF5Vba+qxapaXFhYmNiAkjSrxtkDPgt4dZLbgSuAs5N8bKpTSdIcWDPAVfWeqtpcVVuA84HPVdUbpz6ZJM04zwOWpCYH9ZFEVXUtcO1UJpGkOeMesCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDuqFGNKR6rFTDvENoO66e7KDrGH3Om8P4Nt7HlnX7X3fum7t6OYesCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktRkrJcijz6S/kFgN/BYVS1OcyhJmgcH814QL6mqe6c2iSTNGQ9BSFKTcQNcwFVJdiTZtq8VkmxLspRkaXl5eXITStKMGjfAL6yq5wKvAN6R5CdWr1BV26tqsaoWFxYO8a0BJWmOjBXgqrpz9Oc9wGeAM6c5lCTNgzUDnORJSU7Y+zXwcuCGaQ8mSbNunLMgngZ8Jsne9T9eVX8y1akkaQ6sGeCqug340XWYRZLmiqehSVITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSk4N5Q3bpiHXvGRsP6XonPGV9P9zlG6/ds67bA3jtszas6/auemRdN3dUcw9YkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKajBXgJCcm+b0kX0tyc5Ifn/ZgkjTrxn0p8geAP6mq1yR5PPDEKc4kSXNhzQAn+V7gJ4A3A1TVd4HvTncsSZp94xyCeCawDPxOkq8m+XCSJ015LkmaeeME+BjgucBvVtUZwMPAJatXSrItyVKSpeXl5QmP2ev004eLJE3SOMeAdwG7qupLo+9/j30EuKq2A9sBFhcXa2ITHgEuvbR7AkmzaM094Kq6C7gjydbRj14K3DTVqSRpDox7FsQ7gctHZ0DcBrxleiNJ0nwYK8BVtRNY348OkKQZ5yvhJKmJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpSaom/745SZaBr0/8hjUNz6iqhe4hpHk0lQBLktbmIQhJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpyf8De7/H7kLW/IUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, asgd_pi, 'semi-dual : OT matrix ASGD')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Sinkhorn results\n",
    "---------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEc5JREFUeJzt3X2QXQV9xvHnMYRBDII0OxYIuBY1DmMl4IovKKUwYoIW246jUl+Ktc3YWgdaWt86baUztdY6No462AAqAkUpQsdBtGAJQ6kQu5FoCSGWUpHwYjalSFAEEp7+cU/sNi57Tzb33l/23u9nZofde84953fD7HfPnj33XicRAGDwnlI9AACMKgIMAEUIMAAUIcAAUIQAA0ARAgwARQgwsBew/Urbm/qw3TfbvqblumfYvnF3l2HuCDCGQhOIf7f9Y9v32z7X9kHNsk/bfrj5eMz249O+/uoAZovt58y2TpJ/SbJ0jtt/he1v2P6h7Qds/6vtFzfbvSTJKXPZLvqPAGPes322pL+W9MeSDpT0UknPknSt7X2TvDPJoiSLJH1I0hd3fp1kRd3kHbb32YP7Pl3SVZI+IelgSYdJOkfSo72Zrvf25PEOGwKMea0J0DmS3p3ka0keT/I9SW+QNC7pLXPY5om2N9t+j+0ttu+z/au2T7X93eYo8wPT1j/O9k22H2zW/aTtfZtlNzSrfbs54n7jtO2/1/b9kj6787bmPkc2+zi2+fpQ21O2T5xh3OdJUpJLk+xI8kiSa5J8p7nv/zt10ByNv9P2fzTzfsq2n+Tf4W9s32j7wGm3fdT2/9j+L9srpt1+qO0vN3PfYft3pi37oO3LbV9s+yFJZzS3XWb787a32d5ge2I3/1fNewQY893LJe0n6YrpNyZ5WNLVkl41x+3+fLPdwyT9maTz1In5iyS9UtKf2n52s+4OSX8gabGkl0k6WdLvNXOc0KxzdHPE/cVp2z9YnSP1lbvM/p+S3ivpYtv7S/qspAuTXD/DnN+VtMP2hbZX2H5Gi8f2WkkvlvRCdX5QvXr6QttPsX1es/yUJD9sFr1E0qbmcX5E0gXT4v0FSZslHSrp9ZI+ZPukaZt9naTLJR0k6ZLmttOa+x0k6cuSPtli9qFCgDHfLZa0Ncn2GZbd1yyfi8cl/WWSx9WJxGJJH0+yLckGSbdJOlqSkqxLcnOS7c3R999J+qUu239C0p8neTTJI7suTHKepDskrZV0iKQ/mWkjSR6S9ApJUeeHxFRzJPrMWfb94SQPJvm+pDWSlk1btlDSper8cPiVJD+etuyuJOcl2SHpwmauZ9o+XNLxkt6b5CdJ1ks6X9Lbpt33piT/mOSJaY/3xiRXN9u7SM2/5yghwJjvtkpa/CTnFQ9pls/FfzdhkKSdwfjBtOWPSFokSbafZ/uq5o9/D6lznrlb+KeS/KTLOudJeoGkTyR50nO6STYmOSPJkmb9QyWtmmW790/7/Mc7H0fjOeocrZ6T5LEnu9+0MC9q9vdAkm3T1r1Lnd8edrq7xRz7jdr5YQKM+e4mdf7g9OvTb7S9SNIKSf88gBnOlXS7pOcmebqkD0ia8bzqNLO+DGEz/ypJF0j6oO2D2wyS5HZJn1MnxHOxUdLbJX3VdturMu6VdLDtA6bddoSke6aPNsd5hhoBxrzWnJ88R9InbC+3vdD2uKTL1DknedEAxjhA0kOSHrb9fEm/u8vyH0j6hd3c5sclTSb5bUlfkfTpmVay/XzbZ9te0nx9uKTTJd28m/v7qSSXqvND5Ou2j2yx/t2SviHpr2zvZ/uFkt4h6eK5zjAqCDDmvSQfUScYH1UnhGvV+ZX35Nl+de+hP5L0G5K2qXPa4Iu7LP+gpAubqw7e0G1jtl8nabn+L+R/KOlY22+eYfVt6vxxbK3tH6kT3lslnT2Hx/FTSS6U9BeSrmt+oHVzujpXndwr6Up1zm9/fU9mGAXmBdkBoAZHwABQhAADQBECDABFCDAAFBmpi57xsxYvXpzx8fHqMYChsm7duq1JxrqtR4BH3Pj4uCYnJ6vHAIaK7bvarMcpCAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgyD7VA6DYpk3SiSdWT4FRtmyZtGpV9RQlOAIGgCIcAY+6pUul66+vngIYSRwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFHGS6hlQyPY2SZuq5+ixxZK2Vg/RBzyu+WNpkgO6rbTPICbBXm1TkonqIXrJ9uSwPSaJxzWf2J5ssx6nIACgCAEGgCIEGKurB+iDYXxMEo9rPmn1mPgjHAAU4QgYAIoQYAAoQoBHlO3ltjfZvsP2+6rn6QXbn7G9xfat1bP0iu3Dba+xfZvtDbbPrJ6pF2zvZ/ubtr/dPK5zqmfqFdsLbN9i+6pu6xLgEWR7gaRPSVoh6ShJp9s+qnaqnvicpOXVQ/TYdklnJzlK0kslvWtI/l89KumkJEdLWiZpue2XFs/UK2dK2thmRQI8mo6TdEeSO5M8JukLkl5XPNMeS3KDpAeq5+ilJPcl+Vbz+TZ1vrEPq51qz6Xj4ebLhc3HvL8iwPYSSa+RdH6b9QnwaDpM0t3Tvt6sIfimHna2xyUdI2lt7SS90fyqvl7SFknXJhmGx7VK0nskPdFmZQIMzAO2F0n6kqSzkjxUPU8vJNmRZJmkJZKOs/2C6pn2hO3XStqSZF3b+xDg0XSPpMOnfb2kuQ17IdsL1YnvJUmuqJ6n15I8KGmN5v/5++MlnWb7e+qc1jvJ9sWz3YEAj6Z/k/Rc28+2va+kN0n6cvFMmIFtS7pA0sYkH6uep1dsj9k+qPn8qZJeJen22qn2TJL3J1mSZFyd76nrkrxltvsQ4BGUZLuk35f0T+r8UeeyJBtqp9pzti+VdJOkpbY3235H9Uw9cLykt6pzNLW++Ti1eqgeOETSGtvfUeeA4NokXS/bGjY8FRkAinAEDABF+vKC7IsXL874+Hg/No0eW7du3dYkY9Vz7KkTT/nwnH6Ve/Xf3tDrUWa15m3HDXR/kpRbBnt26don/sED3eE81pcAj4+Pa3Ky1QvCo5jtu6pnAEYVpyAAoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaBIqwAP4xs4AkC1rgEe4jdwBIBSbY6Ah/INHHfHWWd1PgCgl9q8GM9Mb+D4kl1Xsr1S0kpJOuKII3oy3N5i/frqCQAMo579ES7J6iQTSSbGxub9qxsCQN+1CTBv4AgAfdAmwLyBIwD0QddzwEm22975Bo4LJH1mGN7AEQCqtXpHjCRXS7q6z7MAwEjhmXAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFGn1RAxgb3fd5y+Y0/1OPeHXejzJ7HLn7QPdnyQt4MWx9locAQNAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFugbY9mdsb7F96yAGAoBR0eYI+HOSlvd5DgAYOV0DnOQGSQ8MYBYAGCmcAwaAIj0LsO2VtidtT05NTfVqswAwtHoW4CSrk0wkmRjj9UcBoCtOQQBAkTaXoV0q6SZJS21vtv2O/o8FAMOv61sSJTl9EIMAwKjhFAQAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABTp+kw4YD5Y8ZyXz+l+379o/x5PMrtH7p0Y6P4k6bnvXjvwfaIdjoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIm3eFflw22ts32Z7g+0zBzEYAAy7Nq8FsV3S2Um+ZfsASetsX5vktj7PBgBDresRcJL7knyr+XybpI2SDuv3YAAw7HbrHLDtcUnHSPqZl1eyvdL2pO3Jqamp3kwHAEOsdYBtL5L0JUlnJXlo1+VJVieZSDIxNjbWyxkBYCi1CrDtherE95IkV/R3JAAYDW2ugrCkCyRtTPKx/o8EAKOhzRHw8ZLeKukk2+ubj1P7PBcADL2ul6EluVGSBzALAIwUngkHAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQJE2rwcM7PV+8sqj5nS/Ay8b7LfAz/3WDwa6P+zdOAIGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAirR5V+T9bH/T9rdtb7B9ziAGA4Bh1+aJ8I9KOinJw7YXSrrR9leT3Nzn2QBgqLV5V+RIerj5cmHzkX4OBQCjoNU5YNsLbK+XtEXStUnWzrDOStuTtienpqZ6PScADJ1WAU6yI8kySUskHWf7BTOsszrJRJKJsbGxXs8JAENnt66CSPKgpDWSlvdnHAAYHW2ughizfVDz+VMlvUrS7f0eDACGXZurIA6RdKHtBeoE+7IkV/V3LAAYfm2ugviOpGMGMAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAkTbPhAP2evtvuG9O99v3nnt7PMns9rl5yUD3J0lfuXf9wPeJdjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIq0DrDtBbZvsc0bcgJAD+zOEfCZkjb2axAAGDWtAmx7iaTXSDq/v+MAwOhoewS8StJ7JD3xZCvYXml70vbk1NRUT4YDgGHWNcC2XytpS5J1s62XZHWSiSQTY2NjPRsQAIZVmyPg4yWdZvt7kr4g6STbF/d1KgAYAV0DnOT9SZYkGZf0JknXJXlL3ycDgCHHdcAAUGS33pIoyfWSru/LJAAwYjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaDIbj0RA9hbPf6sub0AlO+5t8eTzG7HPfcNdH+S9MMnHhno/p4x0L3NbxwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEVaPRW5eUv6bZJ2SNqeZKKfQwHAKNid14L45SRb+zYJAIwYTkEAQJG2AY6ka2yvs71yphVsr7Q9aXtyamqqdxMCwJBqG+BXJDlW0gpJ77J9wq4rJFmdZCLJxNjY3F4aEABGSasAJ7mn+e8WSVdKOq6fQwHAKOgaYNtPs33Azs8lnSLp1n4PBgDDrs1VEM+UdKXtnev/fZKv9XUqABgBXQOc5E5JRw9gFgAYKVyGBgBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARXbnBdmBvdbWX3zqnO534KIX9XiS2d39m9sHuj9JeuORCwa6v2seGeju5jWOgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoEirANs+yPbltm+3vdH2y/o9GAAMu7ZPRf64pK8leb3tfSXt38eZAGAkdA2w7QMlnSDpDElK8pikx/o7FgAMvzanIJ4taUrSZ23fYvt820/r81wAMPTaBHgfScdKOjfJMZJ+JOl9u65ke6XtSduTU1NTPR6z1rJlnQ8A6KU254A3S9qcZG3z9eWaIcBJVktaLUkTExPp2YR7gVWrqicAMIy6HgEnuV/S3baXNjedLOm2vk4FACOg7VUQ75Z0SXMFxJ2S3t6/kQBgNLQKcJL1kib6PAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEWc9P51c2xPSbqr5xtGPzwryVj1EMAo6kuAAQDdcQoCAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKPK/bk07WnJikdoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, sinkhorn_pi, 'OT matrix Sinkhorn')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "COMPUTE TRANSPORTATION MATRIX FOR DUAL PROBLEM\n",
    "############################################################################\n",
    "############################################################################\n",
    " SEMICONTINOUS CASE:\n",
    "\n",
    " Sample one general measure a, one discrete measures b for the semicontinous\n",
    " case\n",
    " ---------------------------------------------\n",
    "\n",
    " Define one general measure a, one discrete measures b, the points where\n",
    " are defined the source and the target measures and finally the cost matrix c.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_source = 7\n",
    "n_target = 4\n",
    "reg = 1\n",
    "numItermax = 100000\n",
    "lr = 0.1\n",
    "batch_size = 3\n",
    "log = True\n",
    "\n",
    "a = ot.utils.unif(n_source)\n",
    "b = ot.utils.unif(n_target)\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "X_source = rng.randn(n_source, 2)\n",
    "Y_target = rng.randn(n_target, 2)\n",
    "M = ot.dist(X_source, Y_target)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call the \"SGD\" dual method to find the transportation matrix in the\n",
    "semicontinous case\n",
    "---------------------------------------------\n",
    "\n",
    "Call ot.solve_dual_entropic and plot the results.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.92524245 2.75994495 1.08144666 0.02747421 0.60913832 1.8156535\n",
      " 0.11738177] [0.33905828 0.46705197 1.56941919 4.96075241]\n",
      "[[2.20327995e-02 9.26244184e-02 1.09321230e-02 9.71212784e-08]\n",
      " [1.56579562e-02 1.73985799e-03 1.20373178e-04 2.48153271e-05]\n",
      " [3.49227454e-03 8.05110304e-02 4.44694627e-02 3.42874458e-09]\n",
      " [3.15181548e-02 4.34346087e-02 7.17227024e-02 1.28326090e-05]\n",
      " [6.79336320e-02 5.59136813e-03 5.35899879e-02 2.18675752e-02]\n",
      " [8.02083959e-02 3.60364770e-03 4.97032746e-03 1.14377502e-02]\n",
      " [4.87374362e-02 3.36433325e-02 6.09190548e-02 7.33833971e-05]]\n"
     ]
    }
   ],
   "source": [
    "sgd_dual_pi, log_sgd = ot.stochastic.solve_dual_entropic(a, b, M, reg,\n",
    "                                                         batch_size, numItermax,\n",
    "                                                         lr, log=log)\n",
    "print(log_sgd['alpha'], log_sgd['beta'])\n",
    "print(sgd_dual_pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the results with the Sinkhorn algorithm\n",
    "---------------------------------------------\n",
    "\n",
    "Call the Sinkhorn algorithm from POT\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.55535622e-02 9.96413843e-02 1.76578860e-02 4.31043335e-06]\n",
      " [1.21640742e-01 1.25369034e-02 1.30234529e-03 7.37715259e-03]\n",
      " [3.56096458e-03 7.61460101e-02 6.31500344e-02 1.33788624e-07]\n",
      " [2.61499607e-02 3.34255577e-02 8.28741973e-02 4.07427179e-04]\n",
      " [9.85698720e-03 7.52505948e-04 1.08291770e-02 1.21418473e-01]\n",
      " [2.16947591e-02 9.04086158e-04 1.87228707e-03 1.18386011e-01]\n",
      " [4.15442692e-02 2.65998963e-02 7.23192701e-02 2.39370724e-03]]\n"
     ]
    }
   ],
   "source": [
    "sinkhorn_pi = ot.sinkhorn(a, b, M, reg)\n",
    "print(sinkhorn_pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot  SGD results\n",
    "-----------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEfpJREFUeJzt3X2QXQV9xvHnaYiigjA2q8UkuDg6UMQxOCvGIhZBbXgZrLa1vmALY5uxaIXWqfWlKtp2OrWjTetLbYqKLSgiLx3HUSsI1EIR3ISAQAjY+EIQzUZKTagCgad/3JOZ7ZrsnmzO3d/uvd/PzE723nPOPb+7mXz35Nw3JxEAYO79QvUAADCsCDAAFCHAAFCEAANAEQIMAEUIMAAUIcBY0Gyfb/svOridUduxvV8Xcy1Etl9n+6vVcwwTAgx0wD1/Yvsu2z+1/X3bf2X7sc3yL9ve0Xw9bPuhSZc/3ufZWv1ySXJhkpfNch8vt73B9k9sb7N9le3DJi1/pu2LbE8069xl+8O2lzXLj7f96KSfyRbbF9t+3mzmWSgIMNCNv5e0WtLvSDpQ0kmSTpR0sSQlOSnJAUkOkHShpA/supzkjVVD77IvR/62nyHpnyW9VdJBkg6T9FFJj0xafoOkH0g6OskTJR0r6b8kvXDSTf2g+fkcKGmlpDsk/YftE2c723xHgLGg2D7a9nrb221/TtL+k5adYfvaKeunCYBsn2L7puYI7G7b53Y00zMlnSXpdUmuT7IzyW2SfkPSKtsnzOI2z7B9ne2/tX2/7c22f6W5/m7bW23/7qT1p7tvX2/+vL85unzBlNv/saRzJ//8mn1ts728ufwc2/9t+4jdjLtC0neSfC0925NcmuT7zfJzJV2X5I+TbJGkJFuTrEly0dQba25jS5L3SDpP0l/v7c9voSDAWDBsP0bSv0r6F0lPkvR59SLX1gPqHaEeLOkUSX9g+9db7vtjtj+2h8UnStqS5MbJVya5W9I3JL10L2ac7PmSbpH0i5I+I+kiSc+T9AxJp0v6iO0DmnWnu28vav48uDnivn7S7W+W9BRJfzll9v+U9I+SPm37cZIukPTuJHfsZs71ko5oYv7iSTPt8hJJl+71ve+5TNJzbT9hltvPawQYC8lKSYslrUnycJJLJH2z7cZJrknyrSSPJrlF0mcl/WrLbc9KctYeFi+RdO8elt3bLJ+N7yT5VJJHJH1O0nJJ70/yYJKvSnpIvRjP9r79IMmHmyP2n+5m+bnqnVK4UdI96p1W+DlJNks6XtJS9U65bGseHN0V4iWSfrhrfdtvbo7qd9j+p5lmlGT1frEMHAKMheSpku7J/38Hqe+13dj2821f3TwQ9D+S3qjZx3GybZIO2cOyQ5rls/GjSd//VJKSTL3uAGnW9+3u6RYmeVjS+ZKOkvTBKT/3qet+I8mrkoxIOk69o+53NYt/rEk/nyQfSXKwpDXq/UKdzlJJkXT/DOstSAQYC8m9kpba9qTrDp30/QOSHr/rgu1fmrL9ZyR9QdLyJAdJ+rh6R1f76ipJy20fM/nK5vzpSklf62AfM5nuvu0pnNO+FaLtpZLeK+lTkj646xkdM0nyTfVOHRzVXPU1Sa9ss+1uvELS+iQPzHL7eY0AYyG5XtJOSW+xvdj2KyVNjt7Nkp5le4Xt/dX7L/RkB0q6L8nPmli+touhktypXvAutL3S9iLbz1LvvOeVSa7sYj8zmO6+TUh6VNLT295Y80vufEmfkPQG9X75/fke1n2h7d+3/eTm8hGSTlPv/LfU+3s4zvaHmqjL9hJJv7ynfdteavu9kn5P0jvbzr3QEGAsGEkeUu9I6gxJ90n6bfWOtHYtv1PS+yVdKekuSddOuYmzJL3f9nZJ71HzFLE2bH/c0z9f983qPWJ/gaQdkr4i6Rrt3YOE+2KP9y3J/6r3INt1zbnXlS1u7y2SnqzeA2+RdKakM20ft5t171cvuN+yveu+Xy7pA83+71TvAb9lkm5uZrxOvfO77550O09ttt+h3rn9Z0s6vjnfPZDMG7IDQA2OgAGgCAEGgCIEGACKEGAAKDK0b72HniVLlmR0dLR6DGCgrFu3blvzopRpEeAhNzo6qvHx8eoxgIFiu9UrNDkFAQBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAU2a96ABTbtEk6/vjqKTDMVqyQ1qypnqIER8AAUIQj4GF3+OHSNddUTwEMJY6AAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAijhJ9QwoZHu7pE3Vc3RsiaRt1UP0Afdr4Tg8yYEzrbTfXEyCeW1TkrHqIbpke3zQ7pPE/VpIbI+3WY9TEABQhAADQBECjLXVA/TBIN4nifu1kLS6TzwIBwBFOAIGgCIEGACKEOAhZXuV7U22v2377dXzdMH2J21vtX1r9Sxdsb3c9tW2b7d9m+2zq2fqgu39bd9o++bmfr2veqau2F5k+ybbX5xpXQI8hGwvkvRRSSdJOlLSa2wfWTtVJ86XtKp6iI7tlPTWJEdKWinpTQPyd/WgpBOSPEfSCkmrbK8snqkrZ0va2GZFAjycjpH07SSbkzwk6SJJLy+eaZ8l+bqk+6rn6FKSe5Osb77frt4/7KW1U+279OxoLi5uvhb8MwJsL5N0iqTz2qxPgIfTUkl3T7q8RQPwj3rQ2R6VdLSkG2on6UbzX/UNkrZKuiLJINyvNZLeJunRNisTYGABsH2ApEslnZPkJ9XzdCHJI0lWSFom6RjbR1XPtC9snyppa5J1bbchwMPpHknLJ11e1lyHecj2YvXie2GSy6rn6VqS+yVdrYV//v5YSafZ/q56p/VOsH3BdBsQ4OH0TUnPtH2Y7cdIerWkLxTPhN2wbUmfkLQxyYeq5+mK7RHbBzffP07SSyXdUTvVvknyjiTLkoyq92/qqiSnT7cNAR5CSXZKerOkf1PvQZ2Lk9xWO9W+s/1ZSddLOtz2FttvqJ6pA8dKer16R1Mbmq+Tq4fqwCGSrrZ9i3oHBFckmfFpW4OGlyIDQBGOgAGgSF/ekH3JkiUZHR3tx02jY+vWrduWZKR6jn113Gl/M6v/yv372rl9I66TT/ytOd2fJD2y8a453d8Vj37ec7rDBawvAR4dHdX4eKs3hEcx29+rngEYVpyCAIAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIq0CPIgf4AgA1WYM8AB/gCMAlGpzBDyQH+C4N845p/cFAF1q82Y8u/sAx+dPXcn2akmrJenQQw/tZLj5YsOG6gkADKLOHoRLsjbJWJKxkZEF/+6GANB3bQLMBzgCQB+0CTAf4AgAfTDjOeAkO23v+gDHRZI+OQgf4AgA1Vp9IkaSL0n6Up9nAYChwivhAKAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgACjS6oUYwHz32B8/OKvtVj3tmI4nmV4evmtO94f5jSNgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoMiMAbb9Sdtbbd86FwMBwLBocwR8vqRVfZ4DAIbOjAFO8nVJ983BLAAwVDgHDABFOguw7dW2x22PT0xMdHWzADCwOgtwkrVJxpKMjYyMdHWzADCwOAUBAEXaPA3ts5Kul3S47S2239D/sQBg8M34kURJXjMXgwDAsOEUBAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFJnxlXDAQrDo1s2z2u7M2zZ1PMn03nX5a+d0f5L09D+9fs73iXY4AgaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKtPlU5OW2r7Z9u+3bbJ89F4MBwKBr814QOyW9Ncl62wdKWmf7iiS393k2ABhoMx4BJ7k3yfrm++2SNkpa2u/BAGDQ7dU5YNujko6WdMNulq22PW57fGJiopvpAGCAtQ6w7QMkXSrpnCQ/mbo8ydokY0nGRkZGupwRAAZSqwDbXqxefC9Mcll/RwKA4dDmWRCW9AlJG5N8qP8jAcBwaHMEfKyk10s6wfaG5uvkPs8FAANvxqehJblWkudgFgAYKrwSDgCKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAibd4PGJj3fvS6o2a13Z9dNrvtZuukl4zP6f4kadOc7xFtcQQMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFGnzqcj7277R9s22b7P9vrkYDAAGXZv3gnhQ0glJdtheLOla219O8o0+zwYAA63NpyJH0o7m4uLmK/0cCgCGQatzwLYX2d4gaaukK5LcsJt1Vtsetz0+MTHR9ZwAMHBaBTjJI0lWSFom6RjbP/cefknWJhlLMjYyMtL1nAAwcPbqWRBJ7pd0taRV/RkHAIZHm2dBjNg+uPn+cZJeKumOfg8GAIOuzbMgDpH0aduL1Av2xUm+2N+xAGDwtXkWxC2Sjp6DWQBgqPBKOAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKtHklHDDvPeP0O2e13Y5XzO0xyBef+uw53Z8kLT+Vf+bzFUfAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQJHWAba9yPZNtvlATgDowN4cAZ8taWO/BgGAYdMqwLaXSTpF0nn9HQcAhkfbI+A1kt4m6dE9rWB7te1x2+MTExOdDAcAg2zGANs+VdLWJOumWy/J2iRjScZGRkY6GxAABlWbI+BjJZ1m+7uSLpJ0gu0L+joVAAyBGQOc5B1JliUZlfRqSVclOb3vkwHAgON5wABQZK8+qyTJNZKu6cskADBkOAIGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoMhevRADmK9eObJ+Vtt9+mdHdTzJ9I74o+/O6f4kSU9ZMvf7RCscAQNAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFWr0UuflI+u2SHpG0M8lYP4cCgGGwN+8F8eIk2/o2CQAMGU5BAECRtgGOpK/aXmd79e5WsL3a9rjt8YmJie4mBIAB1TbAL0zyXEknSXqT7RdNXSHJ2iRjScZGRkY6HRIABlGrACe5p/lzq6TLJR3Tz6EAYBjMGGDbT7B94K7vJb1M0q39HgwABl2bZ0E8RdLltnet/5kkX+nrVAAwBGYMcJLNkp4zB7MAwFDhaWgAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFBkb96QHZi33nnlq2a13RPPWNTxJNNb/Gtz/5kGTzr1zjnfJ9rhCBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIq0CrDtg21fYvsO2xttv6DfgwHAoGv7UuS/k/SVJL9p+zGSHt/HmQBgKMwYYNsHSXqRpDMkKclDkh7q71gAMPjanII4TNKEpE/Zvsn2ebaf0Oe5AGDgtQnwfpKeK+kfkhwt6QFJb5+6ku3Vtsdtj09MTHQ8Zq0VK3pfANClNueAt0jakuSG5vIl2k2Ak6yVtFaSxsbG0tmE88CaNdUTABhEMx4BJ/mhpLttH95cdaKk2/s6FQAMgbbPgvhDSRc2z4DYLOnM/o0EAMOhVYCTbJA01udZAGCo8Eo4AChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAo4qT7982xPSHpe53fMPrhaUlGqocAhlFfAgwAmBmnIACgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAi/wdOeWKxhqQOygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, sgd_dual_pi, 'dual : OT matrix SGD')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Sinkhorn results\n",
    "---------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEc5JREFUeJzt3X2QXQV9xvHnMYRBDII0OxYIuBY1DmMl4IovKKUwYoIW246jUl+Ktc3YWgdaWt86baUztdY6No462AAqAkUpQsdBtGAJQ6kQu5FoCSGWUpHwYjalSFAEEp7+cU/sNi57Tzb33l/23u9nZofde84953fD7HfPnj33XicRAGDwnlI9AACMKgIMAEUIMAAUIcAAUIQAA0ARAgwARQgwsBew/Urbm/qw3TfbvqblumfYvnF3l2HuCDCGQhOIf7f9Y9v32z7X9kHNsk/bfrj5eMz249O+/uoAZovt58y2TpJ/SbJ0jtt/he1v2P6h7Qds/6vtFzfbvSTJKXPZLvqPAGPes322pL+W9MeSDpT0UknPknSt7X2TvDPJoiSLJH1I0hd3fp1kRd3kHbb32YP7Pl3SVZI+IelgSYdJOkfSo72Zrvf25PEOGwKMea0J0DmS3p3ka0keT/I9SW+QNC7pLXPY5om2N9t+j+0ttu+z/au2T7X93eYo8wPT1j/O9k22H2zW/aTtfZtlNzSrfbs54n7jtO2/1/b9kj6787bmPkc2+zi2+fpQ21O2T5xh3OdJUpJLk+xI8kiSa5J8p7nv/zt10ByNv9P2fzTzfsq2n+Tf4W9s32j7wGm3fdT2/9j+L9srpt1+qO0vN3PfYft3pi37oO3LbV9s+yFJZzS3XWb787a32d5ge2I3/1fNewQY893LJe0n6YrpNyZ5WNLVkl41x+3+fLPdwyT9maTz1In5iyS9UtKf2n52s+4OSX8gabGkl0k6WdLvNXOc0KxzdHPE/cVp2z9YnSP1lbvM/p+S3ivpYtv7S/qspAuTXD/DnN+VtMP2hbZX2H5Gi8f2WkkvlvRCdX5QvXr6QttPsX1es/yUJD9sFr1E0qbmcX5E0gXT4v0FSZslHSrp9ZI+ZPukaZt9naTLJR0k6ZLmttOa+x0k6cuSPtli9qFCgDHfLZa0Ncn2GZbd1yyfi8cl/WWSx9WJxGJJH0+yLckGSbdJOlqSkqxLcnOS7c3R999J+qUu239C0p8neTTJI7suTHKepDskrZV0iKQ/mWkjSR6S9ApJUeeHxFRzJPrMWfb94SQPJvm+pDWSlk1btlDSper8cPiVJD+etuyuJOcl2SHpwmauZ9o+XNLxkt6b5CdJ1ks6X9Lbpt33piT/mOSJaY/3xiRXN9u7SM2/5yghwJjvtkpa/CTnFQ9pls/FfzdhkKSdwfjBtOWPSFokSbafZ/uq5o9/D6lznrlb+KeS/KTLOudJeoGkTyR50nO6STYmOSPJkmb9QyWtmmW790/7/Mc7H0fjOeocrZ6T5LEnu9+0MC9q9vdAkm3T1r1Lnd8edrq7xRz7jdr5YQKM+e4mdf7g9OvTb7S9SNIKSf88gBnOlXS7pOcmebqkD0ia8bzqNLO+DGEz/ypJF0j6oO2D2wyS5HZJn1MnxHOxUdLbJX3VdturMu6VdLDtA6bddoSke6aPNsd5hhoBxrzWnJ88R9InbC+3vdD2uKTL1DknedEAxjhA0kOSHrb9fEm/u8vyH0j6hd3c5sclTSb5bUlfkfTpmVay/XzbZ9te0nx9uKTTJd28m/v7qSSXqvND5Ou2j2yx/t2SviHpr2zvZ/uFkt4h6eK5zjAqCDDmvSQfUScYH1UnhGvV+ZX35Nl+de+hP5L0G5K2qXPa4Iu7LP+gpAubqw7e0G1jtl8nabn+L+R/KOlY22+eYfVt6vxxbK3tH6kT3lslnT2Hx/FTSS6U9BeSrmt+oHVzujpXndwr6Up1zm9/fU9mGAXmBdkBoAZHwABQhAADQBECDABFCDAAFBmpi57xsxYvXpzx8fHqMYChsm7duq1JxrqtR4BH3Pj4uCYnJ6vHAIaK7bvarMcpCAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgyD7VA6DYpk3SiSdWT4FRtmyZtGpV9RQlOAIGgCIcAY+6pUul66+vngIYSRwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFHGS6hlQyPY2SZuq5+ixxZK2Vg/RBzyu+WNpkgO6rbTPICbBXm1TkonqIXrJ9uSwPSaJxzWf2J5ssx6nIACgCAEGgCIEGKurB+iDYXxMEo9rPmn1mPgjHAAU4QgYAIoQYAAoQoBHlO3ltjfZvsP2+6rn6QXbn7G9xfat1bP0iu3Dba+xfZvtDbbPrJ6pF2zvZ/ubtr/dPK5zqmfqFdsLbN9i+6pu6xLgEWR7gaRPSVoh6ShJp9s+qnaqnvicpOXVQ/TYdklnJzlK0kslvWtI/l89KumkJEdLWiZpue2XFs/UK2dK2thmRQI8mo6TdEeSO5M8JukLkl5XPNMeS3KDpAeq5+ilJPcl+Vbz+TZ1vrEPq51qz6Xj4ebLhc3HvL8iwPYSSa+RdH6b9QnwaDpM0t3Tvt6sIfimHna2xyUdI2lt7SS90fyqvl7SFknXJhmGx7VK0nskPdFmZQIMzAO2F0n6kqSzkjxUPU8vJNmRZJmkJZKOs/2C6pn2hO3XStqSZF3b+xDg0XSPpMOnfb2kuQ17IdsL1YnvJUmuqJ6n15I8KGmN5v/5++MlnWb7e+qc1jvJ9sWz3YEAj6Z/k/Rc28+2va+kN0n6cvFMmIFtS7pA0sYkH6uep1dsj9k+qPn8qZJeJen22qn2TJL3J1mSZFyd76nrkrxltvsQ4BGUZLuk35f0T+r8UeeyJBtqp9pzti+VdJOkpbY3235H9Uw9cLykt6pzNLW++Ti1eqgeOETSGtvfUeeA4NokXS/bGjY8FRkAinAEDABF+vKC7IsXL874+Hg/No0eW7du3dYkY9Vz7KkTT/nwnH6Ve/Xf3tDrUWa15m3HDXR/kpRbBnt26don/sED3eE81pcAj4+Pa3Ky1QvCo5jtu6pnAEYVpyAAoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaBIqwAP4xs4AkC1rgEe4jdwBIBSbY6Ah/INHHfHWWd1PgCgl9q8GM9Mb+D4kl1Xsr1S0kpJOuKII3oy3N5i/frqCQAMo579ES7J6iQTSSbGxub9qxsCQN+1CTBv4AgAfdAmwLyBIwD0QddzwEm22975Bo4LJH1mGN7AEQCqtXpHjCRXS7q6z7MAwEjhmXAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFGn1RAxgb3fd5y+Y0/1OPeHXejzJ7HLn7QPdnyQt4MWx9locAQNAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFugbY9mdsb7F96yAGAoBR0eYI+HOSlvd5DgAYOV0DnOQGSQ8MYBYAGCmcAwaAIj0LsO2VtidtT05NTfVqswAwtHoW4CSrk0wkmRjj9UcBoCtOQQBAkTaXoV0q6SZJS21vtv2O/o8FAMOv61sSJTl9EIMAwKjhFAQAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABTp+kw4YD5Y8ZyXz+l+379o/x5PMrtH7p0Y6P4k6bnvXjvwfaIdjoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIm3eFflw22ts32Z7g+0zBzEYAAy7Nq8FsV3S2Um+ZfsASetsX5vktj7PBgBDresRcJL7knyr+XybpI2SDuv3YAAw7HbrHLDtcUnHSPqZl1eyvdL2pO3Jqamp3kwHAEOsdYBtL5L0JUlnJXlo1+VJVieZSDIxNjbWyxkBYCi1CrDtherE95IkV/R3JAAYDW2ugrCkCyRtTPKx/o8EAKOhzRHw8ZLeKukk2+ubj1P7PBcADL2ul6EluVGSBzALAIwUngkHAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQJE2rwcM7PV+8sqj5nS/Ay8b7LfAz/3WDwa6P+zdOAIGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAirR5V+T9bH/T9rdtb7B9ziAGA4Bh1+aJ8I9KOinJw7YXSrrR9leT3Nzn2QBgqLV5V+RIerj5cmHzkX4OBQCjoNU5YNsLbK+XtEXStUnWzrDOStuTtienpqZ6PScADJ1WAU6yI8kySUskHWf7BTOsszrJRJKJsbGxXs8JAENnt66CSPKgpDWSlvdnHAAYHW2ughizfVDz+VMlvUrS7f0eDACGXZurIA6RdKHtBeoE+7IkV/V3LAAYfm2ugviOpGMGMAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAkTbPhAP2evtvuG9O99v3nnt7PMns9rl5yUD3J0lfuXf9wPeJdjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIq0DrDtBbZvsc0bcgJAD+zOEfCZkjb2axAAGDWtAmx7iaTXSDq/v+MAwOhoewS8StJ7JD3xZCvYXml70vbk1NRUT4YDgGHWNcC2XytpS5J1s62XZHWSiSQTY2NjPRsQAIZVmyPg4yWdZvt7kr4g6STbF/d1KgAYAV0DnOT9SZYkGZf0JknXJXlL3ycDgCHHdcAAUGS33pIoyfWSru/LJAAwYjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaDIbj0RA9hbPf6sub0AlO+5t8eTzG7HPfcNdH+S9MMnHhno/p4x0L3NbxwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEVaPRW5eUv6bZJ2SNqeZKKfQwHAKNid14L45SRb+zYJAIwYTkEAQJG2AY6ka2yvs71yphVsr7Q9aXtyamqqdxMCwJBqG+BXJDlW0gpJ77J9wq4rJFmdZCLJxNjY3F4aEABGSasAJ7mn+e8WSVdKOq6fQwHAKOgaYNtPs33Azs8lnSLp1n4PBgDDrs1VEM+UdKXtnev/fZKv9XUqABgBXQOc5E5JRw9gFgAYKVyGBgBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARXbnBdmBvdbWX3zqnO534KIX9XiS2d39m9sHuj9JeuORCwa6v2seGeju5jWOgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoEirANs+yPbltm+3vdH2y/o9GAAMu7ZPRf64pK8leb3tfSXt38eZAGAkdA2w7QMlnSDpDElK8pikx/o7FgAMvzanIJ4taUrSZ23fYvt820/r81wAMPTaBHgfScdKOjfJMZJ+JOl9u65ke6XtSduTU1NTPR6z1rJlnQ8A6KU254A3S9qcZG3z9eWaIcBJVktaLUkTExPp2YR7gVWrqicAMIy6HgEnuV/S3baXNjedLOm2vk4FACOg7VUQ75Z0SXMFxJ2S3t6/kQBgNLQKcJL1kib6PAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEWc9P51c2xPSbqr5xtGPzwryVj1EMAo6kuAAQDdcQoCAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKPK/bk07WnJikdoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, sinkhorn_pi, 'OT matrix Sinkhorn')\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
